!Burgers Equation

program Burgers
    implicit none
    !define numerical precision
    integer, parameter::pre = selected_real_kind(13)
    real, parameter::pi = 3.1415926

    !define cell
    type geo_mesh
        real(pre)::u !state value of the next step
        real(pre)::u0 !state value of the present step
        real(pre)::xc!cell center coordinate
        real(pre)::res!f[i+1/2]-f[i-1/2]
        real(pre)::du!only for reconstruction
        real(pre)::k!only for reconstruction
    end type geo_mesh

    !define local variables
    type(geo_mesh), allocatable ::cell(:)
    real(pre)::xmin, xmax
    real(pre)::dx, dt
    real(pre)::t, tf
    real(pre)::dwl,dwr
    integer::ncells
    integer::nsteps
    integer::itime
    integer::istage
    integer i, j
    !real(pre)::flux
    real(pre)::res
    real(pre)::du, k

    !calculate area
    ncells = 2000
    tf = 5.0_pre
    xmin = 0.0_pre
    xmax = 10.0_pre

    !define size of cell
    allocate (cell(0:ncells + 1))!2 ghost cells
    dx = (xmax - xmin)/ncells

    !initial condition
    do i = 0, ncells + 1
        cell(i)%xc = 0.0025 + (i - 1)*dx
        if (cell(i)%xc >= 0.0 .and. cell(i)%xc <= 0.975) then
            cell(i)%u0 = 1 - cos(2*pi*cell(i)%xc)
        else
            cell(i)%u0 = 0.0
        end if
        cell(i)%u = cell(i)%u0 !1st step
    end do

    !time advancing
    t = 0.0
    nsteps = 0
    time_advancing: do itime = 1, 5000 !general control, max step <= itime
        call output(ncells, itime)
        if (t == tf) exit
        dt = 0.002
        if (t + dt > tf) then
            dt = tf - t
        end if
        t = t + dt
        nsteps = nsteps + 1

! RK2--get du with minmod
        rk_2: do istage = 1, 2
            reconstruction: do j = 1, ncells
                dwl = cell(j)%u - cell(j - 1)%u
                dwr = cell(j + 1)%u - cell(j)%u
                cell(j)%du = minmod(dwl, dwr)
            end do reconstruction

!res initialize
            res_init: do j = 1, ncells
                cell(j)%res = 0.0
            end do res_init

!res computing
            flux_cal: do j = 2, ncells - 1
                cell(j)%res = flux(cell(j - 1)%u, cell(j)%u, cell(j + 1)%u) !inner cells
            end do flux_cal
            cell(1)%res = flux(0.0, cell(i)%u, cell(2)%u)
            cell(ncells)%res = flux(cell(ncells - 1)%u, cell(ncells)%u, 0.0)

!update
            if (istage == 1) then
                step1: do j = 1, ncells
                    cell(j)%u0 = cell(j)%u
                    cell(j)%u = cell(j)%u - (dt/dx)*cell(j)%res
                end do step1
            else
                step2: do
                    cell(j)%u = cell(j)%u - (dt/dx)*cell(j)%res
                    cell(j)%u = 0.5*(cell(j)%u0 + cell(j)%u)
                end do step2
            end if

        end do rk_2
    end do time_advancing

    write (*, *)
    write (*, *) "final time t(sec) = ", t, "by", nsteps, "time steps"
    deallocate (cell)
    stop

contains

    function minmod(a, b)
        implicit none
        real(pre), intent(in)::a, b
        real(pre)::minmod
        if (abs(a) < abs(b)) then
            minmod = a
        else
            minmod = b
        end if
        return
    end function minmod

    function flux(u, v, w)
        implicit none
        real(pre), intent(in)::u, v, w
        real(pre)::flux, k
        real(pre)::dul, dur, du
        dul = abs(v - u)
        dur = abs(w - v)
        du = minmod(dul, dur)
        k = du/dx !不知dx是否默认使用全局变量
        flux = (v + dx*0.5*k)*(v + dx*0.5*k)
        return
    end function flux

    subroutine output(ncells, N)
        implicit none
        integer, intent(in)::ncells
        integer::i, os, N
        integer::j
        open (unit=itime, file="solution.plt", form="formatted", status="unknown", iostat=os)
        write (itime, *) "variables=""x"",""u""/"
        Write (itime, *) "Zone T=""Burgers_RK2 step---""", itime, "/"
        do i = 1, ncells
            write (itime, *) cell(j)%xc, cell(j)%u
        end do
        close (itime)
        return
    end subroutine

end program
